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The normal field instability in magnetic liquids is investigated experimentally by means 
of a radioscopic technique which allows a precise measurement of the surface topography. 

rf\ ■ The dependence of the topography on the magnetic field is compared to results obtained 

^- | by numerical simulations via the finite element method. Quantitative agreement has been 

ON ■ found for the critical field of the instability, the scaling of the pattern amplitude and the 

detailed shape of the magnetic spikes. The fundamental Fourier mode approximates the 
shape to within 10 % accuracy for a range of up to 40 % of the bifurcation parame- 

f^ | ter of this subcritical bifurcation. The measured control parameter dependence of the 

\Q • wavenumber differs qualitatively from analytical predictions obtained by minimization 

^D \ of the free energy. 

•i-H 
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1. Introduction 

• i-H ■ 

Pattern formation has mostly been investigated in syste ms driven far from equilibriu m, 
like Rayleigh-Benar d convection, Taylq r-Couette flow (cf. ICross fc Hohenberall993(l or 



current instabilities IjPeinke et oZJIl992^ . This lopsided orientation is partly due to the 
belief that mainly systems far from equilibrium ca n bring us a step forward to comprise 
fundamental riddles like the origin of life on earth <|Prigoginel ll988'l . 

However, conservative systems may also exhibit the formation of patt erns, a ty pical 
exam ple being given by elastic shells under a buckling load (Taylor 1933: L ange k. Ne well 
Il97l|) . In particular, these n on-dissipative systems have recently gained interest within 
the context of life: following IShinman fc Newell! I|2004h they can describe pattern for- 
mation in plants. Clos ely related ar e surface instabilities of dielectric liquids in elec- 
tric fields IITavlor fc McEwan 1 19651) and its magnetic counterpart, first observed by 
ICowlev fc Rosensweid (J1967J). In comparison with shell structures, the surface instabili- 
ties are experimentally more accessible because the external field can serve as a convenient 
control parameter. 

In order to observe the Rosensweig, or normal field instability, a horizontally extended 
layer of magnetic fluid is placed in a magnetic field oriented normally to the flat fluid 
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surface. When exceeding a critical value B c of the applied magnetic induction, one can 
observe a hysteretic transition between the flat surface and a hexagonal pattern of liquid 
crests. The transition gains some complexity from the fact that under variation of the 
surface three terms of the energy are varying, namely the hydrostatic energy (determined 
via the height variation of the liquid layer), the surface energy and the magnetic field 
energy. As the surface profile deviates from the flat reference state, the first two terms 
grow whereas the magnetic contribution decreases. At the critical field, the overall energy 
is minimized by a static surface pattern o f finite amplitude. 

Usin g the energy minimization principle. lGailitisI l|l969lll977|) and lKuznetsov fc Spektorl 
l|1976|) were able to deduce an amplitude equation that connects the instability to a sub- 
critical bifurcation. However, this result is lim ited to tiny susceptibilitie s xo ^ 1- Only 
recently this deficiency was partly overco me bv lFriedrichs fc Eng cl (2001) for the normal 
field instability, and bv iFriedrichsl l|2002h for the more general case of the tilted field in- 
stability. Besides these achievements, the nonlinear stability of surface p atterns under the 
assum ption of a linear magnetization law was studied numerically by iBoudouvis et all 
dl987l) . 

From an experimental point of view, the above predictions have been poorly inves- 
tigated. This is mainly due to the experimental difficulties in measuring the surface 
profiles. Du e to their co lloidal nature, magnetic fluids are opaque and have a very poor 
reflectivity IjRosensw cia 1985J). They appear black to the naked eye. Thus stan dard op- 
tical techniques such as holography or triangulation IjPerlin. Lin fc Tindll993|) are not 
successful. Moreover, the fully developed crests are much too steep to be measured with 
optical shadowgraphy, as proposed bv lBrowaevs et all Ijl999l) . utilizing the slightly de- 
formed surface as a (de) focusing mirror for a parallel beam of light. Another method, 
recently accomplished bv lWernet et all 1)200 1|) . analyzes the reflections of a narrow laser 
beam in a Faraday experiment. However, it yields onl y the local surface slop e, but not the 
local surface height. This has now been overcome bv lMegalios et all (2005). By adapting 



the focus of a laser beam, a height detection is possible. However, the method has defi- 
ciencies in accuracy, because the beam partly penetrates the magnetic fluid. Moreover, it 
is limited to maxima and minima whereas measurements of the full surface topography 
are not possible. 

A surface profile can be obtained by s i mple la teral observation of the instability, 
as implemented e . g. by iMahr fc R.ehberd lll998rJl for a single Rosensweig peak and 
bv iBacri fc Salinl l|l984h and IMahr fc Rehberel (fl998foh for a chain of peaks, ft has, 
however, the severe disadvantage that the observed crests are next to the container edge. 
Therefore, they are deeply affected by the meniscus and field inhomogeneities. A com- 
parison with the theory for an infinitely extended layer of magnetic fluid has remained 
unsatisfactory for a long time. Consequently, onl y the 'fl at aspects' of the pattern 
like t he wavenumber llAbou. Wesfreid fc , Rouxl 1200 li lLange. Reimann fc Richterl 
2000; Reimann. Richter. Rehber g fc Langd l200.1l) or the dispersion relation 
IIMahr. Groisman fc RchbergJ Il996t iBrowaevs. Bacri. Flament. Neveu fc Perzv nski 
1999J) has been thoroughly investigated in experiments so far. 

Recently we have developed a radioscopic measurement technique that utilizes the 
attenuation of an X-ray beam to measure the full t hree-dimensional surface profile of the 
magnetic fluid far away from the container edges IjR.ichter fc Blasingll200ll) . 

The aim of this article is to compare thoroughly the measured surface profile of the 
instability with the analyt ical predictions and with novel numerical results obtained by 
iMatthies fc Tobiskal <|2005 1 via finite element methods. 

The paper is organized as follows. In the next section we give a description of the 
experimental methods. Then the numerical methods of the corresponding simulations 
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Figure 1 . Setup of the radioscopic measurement of the surface topography 



arc explained. In the fourth section the experimental and numerical results are presented 
in parallel. Eventually, we compare and discuss these results. 



2. Experimental methods 

In the following, we give a description of the experimental setup, the necessary image 
corrections, and the calibration of the height profile. Finally we characterize the utilized 
magnetic liquid. 

We measure the surface topography of the pattern via the attenuation of X-rays passing 
the magnetic fluid layer in vertical direction. Figure Q] shows the experimental setup. A 
container is placed in the centre of a Helmholtz pair of coils. The container is machined 
from Macrolon® and has a diameter of 170 mm, and a depth of 25 mm. In order to 
minimize the field inhomogeneity at the container edge due to the discontinuity of the 
magnetization, we have introduced a 'beach'. The floor of the container is flat within a 
diameter of 130 mm, outside of which it is inclined upwards at 32 degrees, so that the 
thickness of the fluid layer smoothly decreases down to zero towards the side of the vessel. 
The container is filled 10 mm deep with magnetic fluid. This fill ing depth is i n the range 
of the critical wavelength (9.98 mm) which ensures according to lLangel l|200l|) that finite 
size effects in vertical direction can be excluded. 

The two Helmholtz coils have an inner bore of 29 cm (Oswald Magnettechnik) and a 
vertical separation distance of 18.5 cm. The field homogeneity within the empty coils is 
better than 0.5 % within the volume covered by the vessel. The coils are supplied by a 
high precision constant current source (Heinzinger PTNhp 32-40) which allows to control 
the magnetic induction in steps of 1 iiT up to a maximum induction of 40 mT. 

An X-ray tube with a focus of 0.4 mm x 1.2 mm is mounted above the centre of the 
vessel at a distance of 170 cm. The tube has a wolfram anode in order to emit purely 
continuous radiation in the range of the applied acceleration voltage (20 kV to 60 kV) . 
This ensures that the beam hardness can be adapted to the absorption coefficient of the 
fluid investigated. The radiation transmitted through the fluid layer and the bottom of the 
vessel is recorded by means of an X-ray sensitive photodiode array detector. The detector 
provides a dynamic range of 16 bits, a lateral resolution of 0.4 mm, and a maximum frame 
rate o f 7.5 pictures per second. For more details see the article by iRichter fc Biasing 
lJ200lh . 
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Figure 2. Fluid height versus radiation intensity. The inset shows a schematic view of the wedge 
filled with ferrofluid. The relative amount of the X-rays passing through the ramp can be fitted 
nicely with a sum of exponential functions (solid line, equation l2.31 . A simple exponential decay 
with one attenuation coefficient is not sufficient (dashed line, equation 12. 2H . 



2.1. Image improvement 

Because of deficiencies of the X-ray detector, various digital filters had to be applied to 
improve the image quality. 

First, there are so called 'bad pixels'. These are non- functional sensors which must be 
ignored during image processing. The position of those bad pixels has been extracted 
from an empty test image by looking for all pixels the value of which deviates more than 
10 % from the median of their neighbours. There are single isolated bad pixels as well 
as few unusable rows and columns. Correction is achieved by discarding the value of the 
bad pixels and filling the gap with a bilinear interpolation from the neighbourhood. This 
correction of bad pixels can be see n as the first approximation of the Voronoi-Allebach 
algorithm fSa uer fc Alleb ach 1987J). It is a reasonable approach when the non-functional 
pixels do not form large clusters, but are rather scattered. 

Second, the detector exhibits a nonlinear response, and individual pixels are not equally 
sensitive. Let / denote the real incoming intensity and r the response from the pixel, 
both normalized to the range [0 ... 1] . Then we fit the inverse response function of every 
individual pixel with a cubic polynomial 



I(r) = air 3 + a^r 2 + a 3 r + a 4 . 



(2.1) 



The reference data for this model comes from 26 empty test images with different illu- 
mination. 

This method has an additional advantage of equalizing spatial inhomogeneities in the 
illumination caused by the nonuniformity of the intensity of the beam. This non-uniform 
illumination is contained in the image set used for the calibration. Thus we can safely 
assume that after the nonlinear correction the value of every pixel is directly proportional 
to the absorption of X-rays at the corresponding point. The assumption, however, is that 
the illumination itself does not change over time. 
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Figure 3. Three-dimensional rendering of the surface pattern at B = 20.1068 mT. The black 
contour lines are 1 mm apart from each other. 
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Figure 4. Logarithmic Fourier space representation of the surface in figure The black 
contour lines are 6dB apart from each other. 
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2.2. Conversion from intensity values to surface height profiles 

The intensity of the X-rays decreases monotonically with the height of the fluid due to 
absorption. The slope of the surface has no influence on the transmitted intensity, since 
X-rays have indices of refraction very close to 1. 

If the radiation would be monochromatic, then the decay would follow an exponential 
law 

I(x) = Ioe-P*, (2.2) 

where (3 is the attenuation coefficient. However, in order to get a sufficient illumination, 
we have to use a polychromatic X-ray source. Hence, the absorption of the radiation 
depends on the wavelength and the decay cannot be described by the simple exponential 
function (|2.2|) : radiation of higher energy ('hard X-rays') is generally absorbed less than 
radiation of lower energy. In order to get a smooth approximation of the attenuation, we 
fit the intensity decay with an overlay of four exponential functions, as shown in figure[2| 

4 4 

J(a:)=Jo53«ie-* B , $>, = ! (2.3) 

The datapoints were obtained by recording the absorption image of a ramp with known 
shape filled with ferrofluid, see the inset in figure [3 To hold the ferrofluid in position, 
we covered the bottom and the side of the wedge with adhesive tape. The absorption in 
this tape lowers the effective intensity of the X-rays in the fluid by about 1 %, which has 
been taken into account for the calculation of the absorption in the fluid. 

As can be seen from figure [21 the equation (|2.3|) is a satisfactory interpolation method 
that allows us to determine the height of the fluid above every point in the image to 
a resolution of up to 10 |a.m. The absolute accuracy is not as good because of practical 
problems. For example, it is difficult to determine the exact position of the wedge either 
mechanically (because of its sharp end) or from the X-ray image due to the lateral 
resolution of 0.4 mm. This means that the absolute error of the fluid height will be 
around ±0.2 mm. However, when measuring the pattern amplitude defined as a difference 
between the maximum and minimum surface elevation in the unit cell, this systematic 
error cancels. 

After applying the corrections from the previous sections, we finally arrive at the 
surface profile. A three-dimensional reconstruction of one of the recorded profiles is shown 
in figure The corresponding Fourier transform is displayed in figure 0] It should be 
noted that we get the Fourier space representation of the surface elevation using the 
radioscopic method, in contrast to other methods where the transform of a photograph 
is used. This allows for interpretation of the Fourier domain data as amplitudes which 
will be exploited later on in this paper. 

2.3. Tracking one single peak 

Both analytical and numerical predictions are based on the assumption that the pattern of 
peaks is periodic in space. But a perfect periodic lattice is not observed in this experiment. 
For example, one can see a grain boundary between two different orientations at the right 
side of the recorded profile in figure [3J Additionally, as mentioned in the introduction, 
peaks in the centre of the container are least affected by the edge. Thus, we select one 
single peak from there. Because the peak floats slightly under variation of the magnetic 
induction, it is further necessary to track it from image to image. 

To accomplish this, we first determine the position of all peaks with subpixel accuracy 
by fitting a paraboloid to the centre of the peaks, and then construct the Voronoi tesse- 
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Figure 5. X-Ray image and Voronoi diagram constructed from the centre of the peaks 



lation from these positions. The region occupied by the considered peak is then defined 
by its Voronoi cell, i.e. e very point in t his region is closer to the considered peak centre 
than to any other peak (Fortune 1995). We track the movement of the considered peak 
on two consecutive images by finding the Voronoi cell on the later image in which the 
previous position of the peak is situated. Figure [5] visualizes the tesselation. Note that 
the Voronoi diagram also exposes the grain boundary made up of a chain of penta-hepta 
defects. 

The paraboloidal fit also yields the maximum level of this single peak. The minimum 
level in the hexagonal grid is reached at the corners of the unit cell. Therefore, their 
arithmetic mean value is taken as the minimum level. The amplitude of the pattern is 
then defined as the difference between the maximum and minimum level. 

2.4. Properties of the magnetic liquid 

We have selected the commercial magnetic fluid APG512a (Lot F083094CX) from Fer- 
rotec Co. because of its convenient wavenumber and its outstanding long-term stability. 
The critical field of two data series that are 5 months apart differs by not more than 1 %. 
This can be attributed in part to its carrier liquid, a synthetic ester, which is commonly 
used as oil for vacuum pumps. The magnetic fluid is based on magnetite. 

The characteristic properties that have an influence on the normal field instability 
can be found in table ^ The density of the fluid has been measured using a buoyancy 
method and the surface tension coefficient has been determined using a commercial ring 
tensiometer. 

The magnetization law M(H) in the interesting range is shown in figure Because 
the fluid is polydisperse, it cannot be expected that the magnetization law is a true 
Langevin function. However, fitting a Langevin function to the i nitial range of the mag - 
netization curve leads to satisfactory results, as used before bv lBrowaevs et all fl999). 
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Figure 6. Magnetization curve of the ferrofluid APG512a from Ferrotec Co. The solid line is a 

fit with the Langevin equation 12.41 . 



Quantity 

Surface tensionf 

Density 

Viscosity 

Saturation magnetization^ 

Initial susceptibility 



a 
P 

V 

Ml 

Xo 



Value 

30.57 ± 

1236 ± 

120 ± 

14590 ± 

1.172 ± 



Error 

0.1 mNm -1 
5kgm~ 3 
5mPas 
100 Am" 1 
0.005 



f The absolute error of the measurement is unkn own. The error given here 

is taken from the analysis by Harkins & Jordan ( 1930) 

t Note that M s is not equal to the real saturation magnetization, but is 

fitted in a way that equation il'i. II approximates M(H) in the initial range 

up to H - lOkArn -1 . 

Table 1. Material parameters of the magnetic fluid. 



The magnetization law reads therefore 



M{H)=M* S coth(7|H|) 



l\H\ 



H 

W\ 



7 



3xo 

m%' 



(2.4) 



This approximation is, of course, only valid in the initial range up to an internal field of 
about i/i n t = 12kAm _1 . The maximum internal field in the simulation is fully contained 
within this range. 

The onset of the instability can be predicted from these material parameters by the 
linear stability analysis according to lRosensweid (|l985), § 7.1. The critical magnetization 
M c of the fluid layer is given by 



Mj. = 



2 
Mo 



1 



'•() 



(gpv) 



1/2 



(2.5a) 



The Surface Topography of a Magnetic Fluid 9 

where 

r> pi r> 

ra = ^/JI^H/fiQ, fl c = 77 and fi t = —. (2.5b) 

Together with M(H) and the jump condition of the magnetic held at the ground of the 
dish, 

B C = H + M(H), (2.6) 

the critical induction can be determined from these implicit equations. For the values 
from tabled we get B c = 17.22 mT. 

The critical wavenumber k c from this analysis is given by 

k c = M. (2.7) 

V a 

With our parameters, this yields k c — 0.629 mm -1 , which corresponds to a wavelength 
A r = 9.98 mm. 



3. Numerical methods 

Our numerical simulation of the normal field instability is based on the coupled system 
of the Maxwell equations and the Navier-Stokes equations together with the Young- 
Laplace equation which represents the force balance at the unknown free interface. In 
this section we want to describe a numerical algorithm for the simulation of the coupled 
system of partial differential equat ions. Simplified numerica l models have been studied 
already, see lBoudouvis et all l|1987|) and lLavrova et all l|2003|) . In this paper, we focus on 
the computation of the peak shapes employing a realistic magnetization curve in th e form 
of the nonlinear Langevin function. This is in contrast to lBoudouvis et all (|1987[) where 
the stability of hexagonal and square pattern was considered for a linear magnetization 
law. 

We consider a horizontally unbounded and infinitely deep layer of ferrofluid. The 
Maxwell equations for the non-conducting ferrofluid reduce to 



curl H = 0, div B = 0. 

The magnetization inside the fluid is assumed to follow equation l|2.4|l while there is no 
magnetization outside the fluid. The Navier-Stokes equations are given by 

P[^- + {u-V)u) -V-T = -pge z , (3.1) 



dt 



where r is the stress tensor given by 



div u = 0, (3.2) 



r = ri (Vtt + Vu T ) - (p+ —H 2 \ I + H®B. 

Here, / denotes the unit tensor, ® is the tensor product, and p = Phyd + Pm is the sum 
of the hydrostatic pressure j»hyd and the fluid-magnetic pressure 

H 

Pm = Ho / M(h) dh. 



In the static case we are considering here (u = 0), the Navier-Stokes equations reduce 
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to 

Vp = - P ge z + noMVH. (3.3) 

The Young-Laplace equation, which represents the force balance at the interface, reads 
as follows 

[\r\]n = cr/Cn, (3.4) 

where a is the coefficient of the surface tension, K, the sum of the principal curvatures 
and [\t\] the jump of stress tensor. 

Integrating the pressure equation l|3.3|) and inserting the Young-Laplace equation (|3.4|) , 
we obtain the relation 

H 

aJC + pgz = p Q J M(h)dh + ty{M-n) 2 - P() . (3.5) 

o 

Here the two terms on the left-hand side characterize the surface energy and the hydro- 
static energy, while the first two terms on the right-hand side capture the energy due to 
the presence of a magnetic fluid. The pressure po in l|3.5(l is a constant reference pressure. 

For the numerical simulation we consider a bounded domain O x (z b , Zt) which is chosen 
in a way that the hexagon Q contains exactly one peak. Furthermore, the boundaries in z- 
direction are assumed to be far away from the interface. Now, the problem is transformed 
into its dimensionless form by using the strength of the applied field for all magnetic 
quantities and a characteristic length scale / which is usually a fixed multiple of the 
wavelength. The domain obtained in this way will be denoted by fi x (z b , z t ). 

The Maxwell equations in dimensionless form read 

curl H = 0, div B = 0, in fi x (z b , z t ). (3.6) 

The first differential equation in l|3.6|l ensures the existence of a scalar magnctostatic 
potential tp so that H — — Vy. Hence, by exploiting the second differential equation 
of (|3.6|l . we get 

-div(/i(a;, \V(p\)V(p) = in Q x (z b , z t ). (3.7) 

The coefficient function fi(x, H) is given by 



fi(x,H) = i M{H) 



H 



x e Ha, 
x £ £If, 



where £If and f2^ are the subdomains of fix (zj,, z t ) which correspond to the regions inside 
and outside the ferrofluid, respectively. The magnetostatic problem Ij3.7|l is a nonlinear 
uniformly elliptic partial equation. The nonlinearity in 13. 7J) is resolved by a fixed-point 
iteration. The partial differential equation (|3.7|l is equipped with the following boundary 
conditions: (p = —zH^ at the bottom boundary z b , tp = —zH^ at the top boundary z t , 
and dip/dn = at the side boundary. Here, H^ and H^ denote the constant strength 
of the magnetic field outside and inside the ferrofluid in the case of an undisturbed 
interface, respectively. Moreover, H^ can be obtained from H ^ by a single algebraic 
equation. 

In the consideration of the Young-Laplace equation we assume that the interface T is 
the graph of a function \I/ : fi — > R, i.e., 

r= {(x,y,z) el 3 : z = *(x,y), (x,y)en}. 
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( a ) (b) 

Figure 7. Two-dimensional mesh (a) and three-dimensional surface mesh (b). 



Now, the sum of the principal curvatures can be written in terms of Vf. We have 



K, 



-div- 



V* 



y/T+W*\ 



Hence, after dividing by a, the Young-Laplace equation is given by 



-div- 



Vi + |v*| 



A^ = F 



in £7, 



(3.8) 



where F contains all magnetic terms and A is the critical wavenumber k c expressed in 
units of the inverse characteristic length scale l~ l . The differential equation is completed 
by the boundary condition Q'F/Sn = due to symmetry. The Young-Laplace equation is 
a nonlinear elliptic equation which is, however, non-uniformly elliptic. The nonlinearity 
is again resolved by a fixed-point iteration. 

Both the magnetostatic problem (|3.7|) and the Young-Laplace equation l|3.8|) are 
solved approximately by finite element methods. The magnetostatic problem is a three- 
dimensional equation which is discretized by continuous piecewise trilinear functions on 
hexahedra. The Young-Laplace equation is a two-dimensional problem. For its discretiza- 
tion, continuous piecewise bilinear functions on quadrilaterals are used. 

Figure[3shows a mesh for the Young-Laplace equation and a three-dimensional surface 
mesh on the peak. 

We have to solve two large systems of nonlinear algebraic equations which correspond to 
the three-dimensional problem for the magnetostatic potential tp and the two-dimensional 
problem for the function u which describes the unknown free surface. In both cases, fixed- 
point iterations have been applied while the arising linear systems of equations were 
solved by a multi-level algorithm in each step of the iteration. In the three-dimensional 
problem for the magnetostatic potential, a geometric multi-level method has been ap- 
plied based on a family of successively refined three-dimensional hexahedral meshes. 
The Young-Laplace equation is solved on quadrilateral mesh which is the projection of 
the three-dimensional surface mesh onto a plane. The arising two-dimensional problems 
were solved by an algebraic (instead of a geometric) multi-level method. The coupling 
of three-dimensional and two-dimensional problems results in quite difficult data struc- 
tures which are needed for the information transfer between the subproblems. These 
data struc tures had to be develo ped and were implemented in the program package 
MooNMD l|John fc Matthiesll2004|) . 

All iteration processes are illustrated in the flow chart shown in figure 03 
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Figure 8. Flow chart. 



Level 

d.o.f. (Young-Laplace equation, 
d.o.f. (magnetic potential) 

Table 2. Number of unknowns (degrees of freedom = d.o.f.) on different refinement levels. 



2 


3 4 


141 


537 2 097 


6 909 


52 089 404 721 



Table shows that the greatest computational costs are caused by the solution of 
the magnetostatic problem (|3.7[) since the three-dimensional problems have many more 
unknowns than the associated two-dimensional problems. Note that the number of un- 
knowns for the three-dimensional problem increases by a factor of 8 in one refinement 
step while the number of unknowns for the two-dimensional problem increases only by a 
factor of 4. 

We carried out numerical simulations with the material parameters given in table ^ 
Figure El shows the peak height as a function of the applied field for different refinement 
levels of the mesh for the underlying finite element method. The given peak height is the 
difference between the highest point on the surface, i.e. at the mid-point of a hexagonal 
cell, and the lowest point, i.e. at one of the corners of the hexagonal cell. Furthermore, 
the theoretical value for the onset of the instability is shown in Figure OH We see that 
the qualitative behaviour is reproduced even on very coarse meshes. Obviously, one gets 
higher peaks on finer meshes. Moreover, we obtain numerically on the finest considered 
mesh a value for the critical magnetic induction which is very close to the theoretical 
value. 



4. Results 

In the experiment we recorded 540 surface profiles in total, increasing and decreasing 
the magnetic induction in a quasistatic manner from 16.7 mT to 20.1 mT at maximum 
in steps of 0.015 mT. Below this range, the surface remains flat apart from attraction 
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Figure 9. Peak height depending on the applied field strength for different refinement levels. 
The thin vertical line represents the critical induction of the instability (from the linear theory, 
see §H1}. 



of the fluid to the boundary of the container. Above this range the hexagonal pattern 
transforms into a square array of peaks which is not considered in the present paper. 

In the next paragraphs, we will compare characteristic properties from these profiles 
to their counterpart from the simulated peaks. First we compare the amplitude of the 
pattern. Then we look at the wavenumber, and finally the full shape is examined. 

4.1. Scaling behaviour of the amplitude 

Figure 1101 shows the amplitude of the pattern as a function of the applied magnetic 
induction, as defined in section I^TTTl The triangles denote the experimental values where 
the size of the symbols has been chosen to approximate the statistical error of the data. 
The open circles represent the result from the simulations. The solid line is a fit of the 
experi mental data with the solution of an amplitude equation from iFriedrichs fc Engell 
(2001). The root of the amplitude equation reads 



1 



A(e)k c = — [b 2 (l>re)± x /b*(\ I -.-> 



AeK 



(4.1) 



where k c = 0.629 mm 1 is the critical wavenumber of equation (|2.7(l and e — (B 2 — 
B 2 )/B 2 is the bifurcation parameter. The fit parameters are 

B c = 16.747 ±0.001 mT, (4.2) 

61 =0.0889 ±0.0001, (4.3) 

b 2 = 0.0873 ±0.0003. (4.4) 

Using these parameters, the analytical function describes the measured data very 
well. However, it contains three adjustable parameters. The coefficients computed by 
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Figure 10. Pattern amplitude as a function of the magnetic induction. Triangles pointing 
upwards (downwards) represent the experimental data for an increasing (decreasing) magnetic 
field. The symbol size approximates the statistical error. For clarity only every 5th point is 
plotted. The connected circles show the simulation data. The thin vertical lines represent the 
numerically found transition points from (to) the flat surface. The thick solid line is a fit of the 
experimental data with the solution of an amplitude equation, the dashed line is representing 
the unstable branch of the solution . The inset reproduces the qualitative behaviour expected 
from the analysis bv lGailitisI <ll977f) . 



iFriedrichs fc Enge l (2001) are not applicable with our material parameters, because the 
expansion is only valid up to a susceptibility xa — 1-05. So the analytical theory has no 
predictive power for our fluid. 

Let us now compare the measured amplitude with the results from the simulation 
which does not use a single adjustable parameter. Qualitatively, the numerical and ex- 
perimental data compare very well. Both curves show a small bistable range, which makes 
it necessary to control the magnetic field in tiny steps to resolve the hysteresis. In the ex- 
periment, the width of the hysteresis loop is 0.17 mT whereas the numerical data expose a 
hysteretic range of 0.06 mT. For a higher concentrated fluid (xo = 2.2), we have recently 
obser ved a larger hysteretic range of 1 mT using the same setup IjRichter fc Barashenkovl 
2005). This indicates that for even smaller susceptibilities additional care has to be taken 
to resolve any hysteresis. 

For both experimental and numerical data, the amplitude jumps at the critical point to 
a height of about 1.5 mm and reaches about 5 mm at the highest field. When decreasing 
the field, the surface pattern vanishes at the induction B* < B c with a sudden drop from 
about 1 mm. Despite the principal similiarity of experimental and numerical results, the 
agreement is not convincing: at corresponding amplitudes the experimental data are 
shifted to lower fields. 

The critical induction seen in the simulations Bf m — 17.25 mT, i.e. the induction at 
which the jump occurs when increasing the field, is in accordance with the theoretical 
value from the linear stability analysis B c = 17.22 mT, see S 12.41 From the experimental 
data, the critical field can be extracted by the fit with equation (|4.1(l . which yields 
B c = 16.747 mT. It deviates only by 3% from the theoretical value. The difference 
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Figure 11. The pattern emerging as the magnetic induction is increased. The induction at the 
centre of the vessel is (a) B = 16.45 mT, (b) B = 16.51 mT, (c) B = 16.56 mT, (d) B = 16.62 mT, 
(e) B = 16.68 mT, (f) B = 16.73 mT. 

between these two thresholds does not lie within the statistical error, which indicates 
that it is a systematic deviation: with the errors given in section T2. 41 the uncertainty of 
the theoretical value B c is about 1.1 %. 

The imperfection induced by the edge of the bounded container has a great influence 
over the emerging pattern. Figure ITT1 shows six consecutive X-ray images for increasing 
magnetic induction. Because of the inhomogeneous magnetic induction over the vessel, 
pattern formation starts at the edge and expands towards the centre until the whole 
surface is covered with peaks. 

Nevertheless, the simulations can be reconciled with the experimental findings if we 
take the shift in the critical field into account. To see this, we plot both curves in a unifying 
diagram, see figure IT2l Instead of the magnetic induction B, we plot the amplitude as 
a function of the bifurcation parameter e = (B 2 — B 2 )/B 2 . Here, the actual B c is used 
for each data set. Now, the simulation data (plotted as a solid line only) match nearly 
perfectly the experimental data represented by the symbols. Slight deviations can be 
found near e = 0: for e > 0, the experimental amplitudes for decreasing and increasing 
field differ while the simulation produces identical results which fall somewhere in between 
these two curves. For e* < e < 0, i.e. in the range where the surface is bistable, the 
experiment already shows a small-amplitude surface pattern for increasing induction 
which should not be there in the ideal case. 

The range of bistability observed in the experiment is about twice as wide as the one 
found in the simulation. Since the transition is not equally sharp due to the imperfection, 
it is not easy to determine the exact range where the surface is bistable. 

4.2. The wavenumber modulus 

Figure 1131 shows the experimental wavenumber for increasing and decreasing magnetic 
field which has been determined in the Fourier space with high precision (figure 0J . 
Subpixel accurracy can be achieved due to the fact, that the Fourier space representation 
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Figure 12. Pattern amplitude as a function of the bifurcation parameter. B c has been chosen 
as 16.747 mT for the experimental data (from equation 14.21 1 and 17.22 mT for the numerical 
data (from the linear theory) . The simulation data are shown as a solid line only, whereas the ex- 
perimental data are represented by the upward (downward) triangles for increasing (decreasing) 
magnetic field. 



is a direct transform of the surface topography, not the transform of a flat photograph. 
Thus it includes the correct amplitudes. As the magnetic induction B is increased, the 
wavenumber first increases and then decreases slowly. When B is reduced afterwards, k 
in creases again, but to sligh tly smaller values; hence, k exhibits a small hysteresis loop. 

iFriedrichs fc Engeli (2001) predict that the critical wavenumber of maximal growth k c 
should always be larger than the wavenumber of the resulting nonlinear pattern. How- 
ever, in our experiment k c was found to be smaller than the wavenumber of the pattern, 
with the difference of the two being less than 3 %. Further, the computed wavenum- 
ber decreases monotonically, as the magnetic induction is increased. In our experiment, 
instea d, it has a maxim um near the cri t ical p oint. This also contradicts the findings 
of Bacri & Safin (1984) and lAbou et all 1J20011) . which report a constant wavenumber. 
Note that the wavenumber measured here should not be confused with the wavenumber of 
maximal growth which is predicted by a linear stability analysis. T his latter wavenumber 
shows a linear increase with B as measured bv LLange et all (2000). During the nonlinear 
stabilization of the pattern, the wavenumber r elaxes to some other value that is induced 
by the boundary conditions l|Lange et alJ l20011. In all previous experiments the container 
had straight vertical edges corresponding to hard boundary conditions, which forces the 
wavenumber to be an integer multiple of the reciprocal container diameter. In contrast 
to that we equipped our container with a ramp that should give more freedom to the 
wa venumber. Such a ram p has been studied extensively e.g. in the context of convection 
bv iRehberg et all l|l987J1 . It is obvious from figure H3l that our ramp permits different 
wavenumbers for the same magnetic induction. Nonetheless the boundary seems to pro- 
vide a soft pinnin g effect, which select s a wav enumber that is not necessarily the preferred 
one computed bv IFriedrichs fc Engeli l|200ll) . 

Because of the rather small difference between the experimentally found wavenumber 
and k c , the wavenumber modulus k has been fixed to the critical value k c — 0.629 mm -1 
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Figure 13. Wavenumber of the experimental peak pattern, as determined from Fourier space. 
Triangles pointing downwards designate decreasing magnetic induction, upward triangles stand 
for rising induction. The dotted vertical line represents the critical induction from equation 14. 21 . 
the dashed line marks the critical wavenumber used in the sim ulations. The solid line is the 
preferred wavenumber taken from the theory bv lFriedrichs fc Engell (J2001D . but computed for a 
different set of parameters (particularly, \o — 0.35), the linear extrapolation is dotted. 



1.5 



ccfco 




S 
< 



■a 0.5 - 

i 
< 

o 

-5-4-3-2-1012345 
Position • k 

(a) 

Figure 14. Comparison of experimentally (symbols) and numerically (solid lines) obtained 
peak profiles at e — (a) and e = 0.35 (b) 




(cf. § 12. 4|) in the simulations. In principle, a numerical ab initio estimation of the 
wavenumber of the pattern for each value of the induction is possible. This involves 
calculating the surface pattern for different preset wavelengths. Afterwards the preferred 
wavenumber can be selected by the minimum of the total free energy of the simulated 
profiles. However, the computational cost of this technique was too high at the present 
stage. Since the deviation of the experimental wavenumber from the critical one is only 
about 3 % at maximum, it is expected that the simulated profiles are a near match of 
the experiment. 
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Figure 15. Analyzing the first three harmonics of the pattern. Figure (a) displays a logarithmic 
greyscale image of the Fourier transform at B — 20.1068 mT. The wavevectors of the extracted 
harmonics are denoted by arrows. Figure (b) shows the corresponding amplitudes of the modes. 
Lines (symbols) mark experimental (computed) data, respectively. Solid thick lines and crosses 
represent basic modes, solid thin lines and circles modes of the form 2fei, dashed lines and 
triangles modes of the form fei — k-2. 



4.3. The shape of the peaks 

Let us now consider the shape of the liquid crests. Figurc lTH presents a direct comparison 
of the measured profile of a single peak selected from the centre of the dish with the shape 
of the peak obtained in the simulations. The diagram displays a diagonal cut through 
the unit cell from one corner to the opposite corner at two representative bifurcation 
parameters, e = and e = 0.35. There is no adjustable parameter in this comparison, 
apart from centering the peak and normalizing to the wavelength. While the theoretical 
and experimental results differ near the critical value because of the imperfection, as 
discussed in section |4*7l1 the data show a perfect match at higher amplitudes (e = 0.35). 
For small amplitudes, it is expected that the pattern can be approximated by the 
dominating Fourier mode, i.e. a hexagonal pattern made up of three cosine waves: 



A(x) = — (coskix 



COSk-2X 



cosk^x) + — . 



Here the wavevectors fei, 2,3 enclose an angle of 120°, sum up to zero and have the 
same magnitude k. To give a quantitative measure to what extent the higher modes 
are important, we fit a hexagonal grid with the next two higher harmonics to both the 
experimental as well as the simulation data. The wavevectors are fei — fc2, &2 — fe.3, fe3 — fei 
with the absolute value |\/l5 and 2fci, 2fc2, 2fc3 with the absolute value 2k. These vectors 
arc displayed in figure H&l fa). Figure E|(b) shows the dependence of the components on 
the bifurcation parameter. As expected, the basic Fourier mode is the largest component 
which contributes over 90 % even for the highest observed amplitude. 

Since there is only a small amount of the higher harmonics, the shape of the peaks is 
mostly determined by the basic mode and thus should not vary much over the measured 
range. This can be seen directly if we rescale both the experimental and numerical data 
in a way that all peaks have the same width and height. In figure^] we plot five different 
normalized peaks from the whole range. While the experimental values seem to match 
perfectly due to the noise, the numerical data exhibit a certain tendency to sharper peaks 
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Figure 16. The experimental (a) and numerically obtained (b) normalized profiles of different 
peaks show that the shape changes only slightly. 

for an increasing field. However, this effect is very small, so the assumption of an invariant 
shape is a good approximation within the range e < 0.44. 



5. Discussion and Conclusion 

We have studied experimentally and numerically the normal field instability in a fer- 
rofluid. The qualit ative features of this bifu rcation can be described very well by the 
nonlinear theory of iFriedrichs &: Engell 1)20011) . A quantitative comparison with this the- 
ory is not possible, because their approximations are only valid up to xo ~ 1, which is 
exceeded by the initial susceptibility of our fluid xo — 1.17. 

For a quantitative comparison we had to calculate the surface topography numerically 
by a finite element method. Our computations agreed with the measured amplitude to 
within 1 %, provided that the uncertainties of the material parameters and geometri- 
cal imperfections were taken into account by matching the critical induction B c . This 
indicates that a fluid as complex as magnetic liquids can be indeed described well by 
the set of three basic equations, the Navier-Stokes equation, the Maxwell equation and 
the Young-Laplace equation. It turned out to be essential to include the experimentally 
obtained magnetization curve. 

In an attempt to measure the preferred wavenumber - namely the one minimizing the 
free energy - in the nonlinear regime, we have introduced somewhat softened boundary 
conditions in the form of a ramp, which allows for smooth variation of the wavenumber 
as a function of the magnetic field. It turned out that this ramp stab ilizes wavenumbers 
above the critical value of k c , while the theory bv IFriedrichs & E ngcl (2001) predicts the 
preferred number to be below k c . Whether a ramp can be constructed which selects the 
preferred wavenumber, remains to be investigated in the future. 

Because the radioscopic measurement technique allows us to investigate the full profile 
of the peaks, we are able to quantify the ratio of the fundamental mode to the higher 
harmonics under variation of B. It turned out, that the higher harmonics contribute 
less than 10 % in a range up to e — 0.44. This is an encouraging result for further 
analytical treatment of the problem in terms of amplitude equations for the first few 
modes. Whether the contribution of higher harmonics remains small for fluids with higher 
susceptibility, where already visual inspection reveals sharper peaks, remains a topic of 
further investigations. 
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